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Abstract. We consider the estimation of the regularization parameter for the simultaneous deblurring of multiple 
noisy images via Tikhonov regularization. We approach the problem in three ways. We first reduce the problem to 
a single-image deblurring for which the regularization parameter can be estimated through a classic generalized 
cross-validation (GCV) method. A modification of this function is used for correcting the undersmoothing typical 
of the original technique. With a second method, we minimize an average least-squares fit to the images and 
define a new GCV function. In the last approach, we use the classical GCV on a single higher-dimensional image 
obtained by concatanating all the images into a single vector. With a reliable estimator of the regularization 
parameter, one can fully exploit the excellent computational characteristics typical of direct deblurring methods, 
which, especially for large images, makes them competitive with the more flexible but much slower iterative 
algorithms. The performance of the techniques is analyzed through numerical experiments. We find that under 
the independent homoscedastic and Gaussian assumptions made on the noise, the three approaches provide almost 
identical results with the first single image providing the practical advantage that no new software is required and 
the same image can be used with other deblurring algorithms. 

Key words. Methods: data analysis - Methods: statistical - Techniques: Image processing 



1. Introduction 

An important problem in image processing is the simulta- 
neous deblurring of a set of observed images (of a fixed ob- 
ject), each of which is degraded by a different point spread 
function (PSF). Consider, for example, images obtained 
by the Large Binocular Telescope (LBT). This instrument 
consists of two 8.4m mirrors on a co mmon mount with a 
spacing of 14.4m between the centers IjAneel et al. 1 1998). 
For a given orientation of the telescope, the diffraction- 
limited resolution along the center-to-center baseline is 
equivalent to that of a 22.8m mirror, while the resolu- 
tion along the perpendicular direction is that of a single 
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8m mirror. One way to obtain an image with improved 
and uniform spatial resolution is to simultaneously decon- 
volve the images taken with different orientations of the 
telescope. 

In a recent paper, IVio et al~l l)2004j) presented an ef- 
ficient solution based on a single image that improves 
the performance of iterative algorit hms developed in the 
conte xt of a least-squares approach ijBertero fc Boccacci I 
2000). An important advantage of iterative algorithms is 
the ease with which constraints such as nonnegativity of 
the solution can be implemented. However, this incurs a 
high computational cost (especially in case of slow conver- 
gence). In addition, there is a lack of a reliable stopping 
criteria. Although these issues may not be critical for im- 



2 



R. Vio, P. Ma, W. Zhong, J. Nagy, L. Tenorio, & W. Wamsteker: GCV for multiple image deblurring 



ages of moderate size, they are critical for large images 
(« 10 7 — 10 8 pixels) or studies that require the deblur- 
ring of a large number of images. In this respect, direct 
methods are potentially more useful. For example, with 
a fixed regularization parameter, Tikhonov deblurring re- 
quires only two two-dimensional discrete Fourier trans- 
forms. 

A limitation of any regularization method, and in par- 
ticular of Tikhonov's approach, is the need to get a com- 
putationally efficient and reliable estimate of the regular- 
ization parameter. We present three different methods for 
such a task in the case of multiple spatially invariant PSFs 
and additive Gaussian noise. In Section [3 we formalize 
the problem and propose three solutions in Sect.El Their 
performance is studied through numerical experiments in 
Sect.0] Section [S] closes with final comments. 

2. Formalization of the problem 

The image restoration problem is to find an estimate 
f(x,y) of a fixed two-dimensional object fo(x,y) given 
p observed images gj(x, y) (j — 1, 2, . . . ,p) each of which 
is degraded by a different spatially invariant blurring op- 
erator; that is, 



where each Kj(x,y) is a known space invariant PSF. 

In practice, the observed images are discrete N x N 
arrays of pixels (images are assumed square for simplicity) 
contaminated by noise. Therefore, model JIJ has to be 
recast in the form 

AT-l 

gj(m,n) = Kj(m — m\n — n')fo(m',n')+'Wj(rn,n), 

m',n'— 

(2) 

where Wj(m,n) will be assumed to be Gaussian white 
noise with standard deviation a w = a w . Note however 
that the methods in Sect. may be used when the noise 
is not white but has a correlation structure that is known 
up to a constant factor. 

If the central peak (i.e., support) of each PSF is much 
smaller than the image and the object does not have struc- 
tures near the image boundaries, then the convolution 
product in Eq. J5J can be well approximated by cyclic 
convolution, which leads to 

gj{m,n) = Kj(m,n)f {m,n) +Wj(m,n), (3) 

where the symbol """" denotes the discrete Fourier trans- 
form (DFT). 

In order to estimate f using Tikhonov regularization, 
it is helpful to use matrix-vector notation. Rewrite Eq. (J2J 
as 

9j .'V/. »\ (4) 

where the arrays g-, f and Wj are obtained by column 
stacking gj(m,ri), fo(m,n), and Wj(m,ri), respectively. 



The matrix Aj is block-circulant and defined by cyclic 
co nvolution with the ?'th PSF. 

iBertero & Boccaccil l|2f)f)f)j) estimate / by Tikhonov 
regularization where the estimate f x minimizes an average 
fit to the images subject to a smoothness constraint: 



fx = argmin f £ \\Ajf - 9] \\ 2 + A 2 || Lf\ 
= argmin (\\Af~g\\ 2 + X 2 \\Lf\\ 2 ), 



(5) 



(6) 



where the global imaging matrix and global image are 
defined, respectively, as 



A: 



( A1 \ 




(9i\ 


A 2 




92 




, Q = 




\a p J 




\9 P J 



(7) 



A is the regularization parameter and L is, usually, the 
identity matrix (for energy bound) or a discrete Lap lacian 
operator (for smoothness constraint) l|Wahba Ill990|) . The 
solution f \ satisfies the normal equations 



X 2 L T L ) f x 



(8) 



Kj(x -x',y- y')f {x', y')dx'dy', (1) where 



A = £AjA J> & = J2 A j9j- 0) 

3=1 3=1 

For a fixed value of A, its explicit solution is given by: 



f x (m, n) = ^ — - , 

lC(m, n) + \ 



(10) 



i=i 



where K,(m,n) = Y^=i l^'( m : n )| 2 - 

In the rest of the paper we consider the selection of 
the regularization parameter A for three different methods 
whose normal equations reduce in each case to ijSjl. 

3. Estimating the regularization parameter 

3.1. A first single-image approach 

The Tikhonov regularization JSJ is transformed into a 
single-image problem © through the use of the higher 
dimensional global image q. We now define an alternative 
single-image approach whose dimensionality is the same 
as that of the original images. 

The solution of the normal equations iJSJ depends on 
the original data q only through -d. In fact, in the case of 
homoscedastic uncorrelated Gaussian noise, ■& and q con- 
tain the same information about f n . This can be just ified 
with a sufficiency argument on i9 l|Lehmann fc Casellal 
1998§). It thus seems natural to pose the problem directly 
through i?: Since the expected value of $ is Af , we pro- 
pose to estimate / by regularizing the misfit || Af — ft \\ 2 . 
But one problem with ft is that its covariance matrix 



R. Vio, P. Ma, W. Zhong, J. Nagy, L. Tenorio, & W. Wamsteker: GCV for multiple image deblurring 



3 



is Cov(#) = a 2 ^ A, which is no longer proportional to 
the identity and therefore the classical generalized cross- 
validation (GCV) is not expected to work well. In fact, our 
simulations show that GCV consistently underestimates A 
when used with i9. However, since A is known, one can 
easily correct for this correlation using a square-root of 
the symmetric matrix A: the image 



« = A- 1 ' 2 ® 



(11) 



A 1 11 

has mean A 1 f and covariance matrix 
of„ A~ 1/2 AA~ 1/2 = a 2 w I. Hence, we define the Tikhonov 
functional 



/ A = argminf||A 1 / 2 /-^|| 2 +A 2 || J L/ 



(12) 



It can be easily shown that the normal equations for 
l|12|) are given by JSJ. This result indicates that the func- 
tional JSJ can be expressed in terms of a single-image <i 

A 1/2 

with PSF A ' , whose Fourier domain representation is 



-/ \ •^K* j {m,n)'gj{m,n) 

am,n) — > — ^ . 

~[ /CV2( m ,n) ' 



(13) 



and K}/ 2 {m, n), respectively. Through f(m,n), the solu- 
tion HI U| l for f\(m,n) can be expressed in the form 



f\(m,ri) 



K}/ 2 (m, n)<f(m, n) 
lC(m, n) + X 2 



(14) 



Since the covariance matrix of <r is proportional to the 
identity, we can make use of the classical generaliz ed cross- 
validation (GCV) fun ction to estimate A (e.g., IWahba I 
ll99(]UVio etal~ll2003h : 



GCV W (A) 



V> ' l- r Tr[H^]/N 2 ^' 1 



(15) 



where H M = A 1/2 ( A + X 2 L*L ) 1 A 1/2 , with frequency 
domain representation 



GCV W) (X) 




<f(m, n)\8{m, n)Y 



K(m, n) + X 2 \S(m, n)p 
\S(m, n)\ 2 \ 



/C(m, n) + X 2 \8{m, n)| 5 



(16) 



Here, S(m, n) provides the two-dimensional DFT of matrix 
L. However, despite its theoretical justification and ade- 
quate practical performance, GCV tends to underestimate 
the regularization parameter (too small a A) in about 10% 
of the cases. A corrected GCV can be used to control this 
without sacrificing its generally good performance. The 



corrected weighted GCV is: 



GCV Q W (A) 



V 2 



A 2 <f(m, n)\6(m, n) 



IC(m, n) + X 2 \8(m, n)| 2 

X 2 \S(m, n)\ 2 



I 



N 2 ^/C(m,n) + A 2 |<5(m,n)|^ 



(17) 



where a > 1. It reduces to (|TB|) when a = 1 and yields 
smoother estimates as a increases. 

Empirical studies suggest good values of a for prac- 
tical use in the range of 1.2 and 1.4. Modifications of 
this so rt have been suggest ed by various researchers 
e.g.. iNvchka et. al. I Il998t). However , only recently 
Cummins et al. l200lHKim fc Gu Il2004[) have systematic 
simulation studies been conducted to evaluate the perfor- 
mance of the corrected GCV in the context of statistical 
nonparametric modeling. The optimality of the corrected 
GCV (|17fl can be interpreted in a way similar to that of 
GCV (|16fl . Instead of estimating a relative loss function, 
Eq. (|17fl estimates a weighted combinati on of the average 
squar ed bias and the average variance l|Cummins et al.1 
l20f)lh . It should also be noted that GCV is a function of 
both A and a, but it is difficult to decide whether GCV 
should be minimized over a. Current research efforts are 
still trying to resolve this issue. 

It is important to emphasize that the problem of un- 
derestimation is not unique to GCV, it is common to vir- 
tually all regularization parameter selection techniques. 
For example, t he generalized maximum likelihood (GML) 
estimate of A IjWahba Ill985|) asymptotically tends to de- 
liver slightly smaller estimates than GCV. An attractive 
feature of our simple correction is that it is not restricted 
to GCV, and sheds a fresh light on underestimation cor- 
rections for GML and other approaches. 

The minimization of the GCV function cannot be done 
analytically, some numerical iteration method has to be 
used. In our simulation we use the Newton-Raphson algo- 
rithm. 



3.2. Multiple-image approach 

Although the single-image <; is adequate for applications 
where the most important information is preserved in the 
single image, here we define GCV for functions J5J and 
©, and compare their performances with that obtained 
with the single-image 

We start by defining a multiple-image ordinary cross 
validation (CV) function. One such function can be ob- 
tained by estimating the prediction error of each image 
using a leave-one-out estimate 



CV(A) 



1 p 1 N 



p 1 — ' iV 2 

3=1 i=l 



(18) 



Here, g^-i is the prediction based on the estimate of f 
that uses all the data except the ith entry in g -. Simple 



4 



R. Vio, P. Ma, W. Zhong, J. Nagy, L. Tenorio, & W. Wamsteker: GCV for multiple image deblurring 



algebra leads to an equation for g 3 \-i in terms of g^ (see 
Appendix 1X1) 



9j,- 



9ji 



1 _ ( H j)ii 



(19) 



where, 



9 j = A jf\, // 



A j {A + X 2 L T L)- 1 A 1 j '. (20) 



The multiple-image CV function can then be written as 

p N 2 



CV(A) 



1 



(21) 



The multiple-image GCV is obtained by modifying (|18fl 
to be rotationaly invariant. That is, the value of the GCV 
function should remain unchanged if each g i is rotated. 
There is no unique way to achieve such rotational invari- 
ance but one method is to substitute the average diagonal 
value Tr [Hj]/N 2 for {Hj) u 



(22) 



3=1 



where r.j = — g^. Note that Eq. reduces to the 
classical GCV when p = 1. 

The implementation of Eq. (1221) is done in the fre- 
quency domain. It can be rewritten as 



GCV^>(A) = ^-£^ 



2 P 



TV, 



i=i 



V; 



where 



Kj(m, n)g(m, n) 



gj(m,n) 



mn K(m,n) + \ 2 \6{m,n)\ 2 
with g(m,n) = Xw=i -^f ( m ; n )di( m > n )i an( i 

\Kj(m, n)| 2 



(24) 



Eh 



/C(m, n) + X 2 \S(m, n)\ 2 



(25) 



To avoid underestimating the regularization parame- 
ter, the corrected GCV for (1221) can be derived as 



]/JV*)' 



(26) 



Its form in the frequency domain is given by Eqs. I]23fl- 
JUJl, and 



EN 



a| i^(m, 



/C(m, n) + A 2 |<5(to, n)| 5 



(27) 



3.3. /4 second single-image approach 

An alternative form for the GCV functional (122[l can 
be obtained by applying the standard GCV directly to 
Eq. © 1 based on the single global-image q. In this case 
the rotation invariance is enforced in a higher dimensional 
space of concatenated images which leads to 



GCVi e) (A) 



1 



pN 2 (l - a Tr[H {s) ]/{pN 2 )) 2 

1 gjlkff 
piV 2 (l-a T.^[H ] ]/{pN 2 )f 



(28) 



where = A(A + \ 2 L t L)~ 1 A t . 

Following the usual int erpretation of the denomina- 
tor in the GCV function (Wahba 1990), we see that 
GCV (9) (A) and GCV (e) (A) differ in that the former 
weights each residual error by the noise degrees of freedom 
of each image; images with a larger trace of Hj receive a 
higher weight. In our case the two GCV functions provide 
identical results (see below). 

It can also be shown that even though <; and q lead 
to the same normal equations, their solutions may differ 
because their GCV may give different answers. In fact, 
one can check that 



GCV W (A) -pGCV (e) (A) = 



N 



(29) 



(23) 4. Some numerical experiments 

We present the results of numerical simulations to test 
the reliability of the regularization parameter estimated 
through the GCV functions defined in Sect. [21 

Figures 11161 compare the deblurring of an extended 
object with a sharp outline and a set of point-like ob- 
jects. We have chosen these particular examples because 
their restoration is a difficult problem. Eight images are 
available for each object. In each case the PSF is a bidi- 
mensional Gaussian with dispersion along the major axis 
set to twelve pixels, and to four pixels along the minor 
axis. Their inclinations take equispaced values in the range 
0° — 160°. Gaussian white noise is added with two differ- 
ent levels of noise by setting the standard deviations to 
1% and 10%, respectively, of the maximum value of the 
blurred images. Three deblurring methods have been used: 
(I) single-image Tikhonov based on the functional H12J1 
and the GCvi function (fT7|b (of = 1.4): (II) Tikhonov 
based on the functional (JSJ) and GCvi 9 ' function $Zjofy 
(a = 1.4); (III) a classic iterative Projected Landweber 
Method applied to the single-image i mages for the itera - 
tive methods obtained as explained in I Vio et al. I l)2004tl . 
that can be used to enforce nonnegativity of the solution 
(a maximum number of 2000 iterations has been adopted, 



We thank Prof. M. Bertero for this suggestion. 



R. Vio, P. Ma, W. Zhong, J. Nagy, L. Tenorio, & W. Wamsteker: GCV for multiple image deblurring 



5 



that is sufficient to provide results that do not show fur- 
ther appreciable visual changes). 

There is no need to test the GCV functional l|28[l be- 
cause in our case it is identical to (|26[1 . The reason is that 
the shape of all the selected PSFs is the same, they only 
differ in their orientation, and since orientation informa- 
tion is lost through the trace, the denominator Dj does 
not depend on the image index j. That is, all the images 
share the same effective degrees of freedom. 

Since the Tikhonov technique does not enforce the non- 
negativity constraint, the deblurred images obtained with 
the single-image approach are also shown with their nega- 
tive values set to zero. Finally, for each of the experiments, 
the true "optimal" value (in the sense of the minimization 
of the rms of the true residuals) for the regularization 
parameter provided by the GCV function (fP7|l has been 
determined, with an accuracy of two significant digits, by 
a simple grid search. In all cases, the "true" and the esti- 
mated regularization parameter have been almost identi- 
cal. In particular, the relative difference between the rms 
of the true residuals of the corresponding solutions was 
always less than 0.3%. 

As shown by the figures, the multiple and single-image 
approaches lead to very similar, if not identical, residual 
rms values. 

As expected, Projected Landweber provides the best 
results in terms of visual appearance and residual rms. In 
fact, as is well known, the nonnegativity constraint limits 
the ringing effect around the sharp structures in the image. 
However, in these experiments the difference seems impor- 
tant only for point-like objects. In any case, these results 
are obtained at a very high computational cost (each it- 
eration requires the computation of two two-dimensional 
DFTs). 

Apart from the nonnegativity constraint, the worst 
performance of the Tikhonov method is due to the 
functional L that tends to smooth out structures with 
sharp features. This characteristic makes Tikhonov more 
suited for the restoration of images of diffuse objects. In 
fact, Figs. 17191 show that for diffuse objects, single-image 
Tikh onov provides competitive restoration quality (see 
also IBertero fc Boccacci Ill998|) . Because of the presence 
of features close to the borders of the images, calcula- 
tions for the image <r(m, n) and the deblurring operation 
have been i mplemented with t he windowing procedure as 
explained in IVio et al. I jjjfljj) . This has required the re- 
moval of a border 50 pixels thick from the image and has 
only been done with the single-image Tikhonov approach 
because, given the necessity of windowing each single im- 
age, it would be tortuous to implement the other version 
of the algorithm. 

In order to test the stability of the A estimates, for each 
experiment 100 realizations of the corresponding noise 
process are generated and added to the image. In gen- 
eral, the estimates of this parameter appear quite stable; 
as measured by the rms of the estimated values divided 
by their median, they fall in the range 1% — 2%. Fig. ^| 
shows the histogram of the values of A obtained with the 



single-image Tikhonov approach in the experiments with 
the worst relative variation. 

5. Concluding remarks 

We have considered the estimation of the regularization 
parameter for the simultaneous Tikhonov deblurring of 
multiple noisy images using three types of GCV func- 
tions based on different data transformations. The first 
approach reduces the multiple-image problem to single- 
image deblurring that allows the use of the classical GCV, 
but we use a corrected GCV that reduces its tendency to 
undersmooth the final solution. A second approach is per- 
haps the most natural as it is a straightforward extension 
of a classical single-image approach. It aims to achieve 
a balance between the smoothness of the unknown im- 
age and its average goodness of fit to the observed im- 
ages. The regularization parameter is estimated through 
a multiple-image extension of the traditional GCV func- 
tion. The third approach is based on a higher dimensional 
single image defined by concatenating all the images into 
a single vector and then using the classical GCV. 

Our numerical experiments indicate that the three 
approaches provide very similar results, with the lower- 
dimensional single image providing the practical advan- 
tage that no new software is required and the same image 
can be used with other deblurring algorithms. 

We have assumed that the noise in all the im- 
ages is independent Gaussian and of constant variance. 
Combinations of the methods we propose can be used 
when some of these assumptions are not met. For example, 
if there are k groups of images with different variances, a 
weighted Tikhonov regularization method can be applied 
to a multiple- image deblurring of k single images. In addi- 
ton, the multiple-image approach can be used if the noise 
in the original image is correlated bu t the covar iance ma- 
trix is known up to a constant (e.g., iGu Il2002h . 

The extreme computational efficiency, coupled with 
the stability and reliabilty of the estimated regularization 
parameters, make Tikhonov's method competitive with 
the more flexible, but much slower, iterative algorithms. 
In fact, and especially for very large images, even in situa- 
tions where these algorithms are expected to provide bet- 
ter results (e.g., when the nonnegativity constraint can be 
exploited), the Tikhonov approach is still a valuable re- 
source (see also IBertero fc Boccacci ll200oTl . For example, 
it can be used to obtain initial solutions for iterative al- 
gorithms that may speed up the convergence. Figure 1111 
shows the residual rms as a function of iteration num- 
ber for the Projected Landweber algorithm started at a 
Tikhonov solution and at the uniform image traditionally 
used. The figure shows that the saving in computation 
time can be considerable (here it is more than a factor of 
10). 
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Ordinal Imaoe 



PS F of tte F irst Image 



Mstfl&l !: rm&= UQ74S72 / 7. = Q.0S7122 r*tannegatwe MeUlod : rma= O07SS7 






First Blurred Image 



First Blurred Image + Norr* 



Mettled III rma-ff.ee.5755- ' 







Fig. 1. Original image, blurred version and blurred + Fig. 2. Deblurring of the image in Fig. ^ with the meth- 
noise version of the first image in the set (see text) and ods decribed in the text: (I) single-image Tikhonov l|12f> 
corresponding PSF. The images are 256 x 256 pixels, and an d GCV^ ) $T7$; its version with the negative elements 



the standard deviation of the noise is 10% of the maximum 
value of the blurred image. 



zeroed; (II) Tikhonov © and GCVL 9j (III) a classic 

iterative Projected Landweber Method (convergence after 
61 iterations). Tikhonov has been used with L = I and 
a = 1.4. 
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Nonnegatf/e MetTiod I: rms= OflS^TfiS 



Original Image 



iP^Ptsttte First Image 





Method I!. ".^7 ' a = 03.24577 



Method III: rms= 0.05051 



First Blurred Image 



First Blurred Image*- Noise 





Fig. 3. The same as Fig. but obtained from an im- Fig. 4. Original image, blurred version and blurred + 

age contaminated with noise whose standard deviation is noise version of the first image in the set (see text) and 

1% of the maximum value of the blurred image. For the corresponding PSF. The images are 256 x 256 pixels, and 

Projected Landweber algorithm 2000 iterations have been the standard deviation of the noise is 10% of the maximum 

carried out. value of the blurred image. 
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= :00061507 / k= 0.05729 



Nonnsgalive Method I: 



= 0.00607 1G / A= 0.0059007 



Nonnegative Method I: 



= 0.0061506 t X= 0.16074 



Method III: rms = 0.0059459 



= 00060716 / 1= 0.016435 



Method III: rms = 0.0059076 



Fig. 5. Deblurring of the image in Fig. 0] with the meth- 
ods described in the text: (I) single-image Tikhonov 03 
and GCV£° 0; its version with the negative elements 
zeroed; (II) Tikhonov © and GCvl 9) (III) a classic 

iterative Projected Landweber Method (2000 iterations). 
Tikhonov has been used with L = I and a = 1.4. 



Fig. 6. The same as Fig. but obtained from an image 
contaminated with noise whose standard deviation is 1% 
of the maximum value of the blurred image. 
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Original Image P£F of trie First Image 




Fig. 7. Original image, blurred version and blurred + 
noise version of the first image in the set (see text) and 
corresponding PSF. The original images was 600 x 600 
pixels, but only the central 500 x 500 pixels are displayed 
(see text). The standard deviation of the noise is 10% of 
the standard deviation of the blurred image. 



irracpc/m.n) Method I: rms = 0.03732 / X = 




hbnnegative Msttod I: rma^O.03729 hfetrndlll: rma= 0.03625 




Fig. 8. Image <r(m, n) for the image shown in Fig. 01 single- 
image Tikhonov deblurring (L = discrete Laplacian), its 
version with the negative elements zeroed, and the iter- 
ative Projected Landweber deblurring (convergence after 
24 iterations). 



10 



R. Vio, P. Ma, W. Zhong, J. Nagy, L. Tenorio, & W. Wamsteker: GCV for multiple image deblurring 



Image ym.n) 





hbnnegative Msttod I: rma^O. 03307 



Msttod III: rms=0.05!S91 





Uniform 

Tikhonov 




200 250 300 
Number of iterations 



Fig. 11. Convergence rates for the Projected Landweber 
algorithm in the numerical experiment corresponding to 
Fig. |21 Two starting guesses are used: the uniform image 
typical of the classic implementation of the algorithm, and 
the nonnegative single-image Tikhonov shown in Fig. |3| 



Fig. 9. The same as Fig. |H1 but obtained from an image 
contaminated with noise whose standard deviation is 1% 
of the standard deviation of the blurred image. Tikhonov 
has been used with L = I. For the Projected Landweber 
algorithm 2000 iterations have been carried out. 

35 I 1 1 1 1 1 1 

30 - 
25 - 
| 20 - 

I 

m 

JD 

1 15- 
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-I 
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1 X10- 3 

Fig. 10. Histogram of the values of the estimated param- 
eter A for the 100 experiments based on Fig. 
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Appendix A: Derivation of Ea. lfTT?)) 

Let ei be the vector with 1 in the ith entry and zero 
everwhere else, and f k _ f the estimate of / that does not 
use the ith entry of g k . We can write 

9k,i = efg k) g k -t = ejA k f k _, t . 

Let Pi be the projection that deletes the ith entry: 

9k,-i = p i9k 
so that Pj Pi = I — e^ef. Let 



a ki = A^Bi, 
then, by definition, 



H k = A k (A + XI)- 1 Al, 



-l 



' I J2 A l9i+AlPjP ig> 

i^k 

= (a + AJ- A T k e ie jA k ^ ^-Ale i ejg k 

_ x (A + \I)-ig k *al(A + \I)-i 



{A + XI) 
It follows that 



1- A + XI)- 1 aki 



del 9 k 
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Akf k ,-i = A k f ~ H k e t efg k 

HkejO^jA + XI yh!> - H k eiefH k eiefg k 
1 - (H k ) u 



T " (-""&) 

e t A kfk,-r - Sfc* = 



' * Kl l-(H k ) u 

x iaKA + Xir^-efg, 

(Hk)i 



1 - (H k )i 

(Hk)u 
1 - [H k )i 



(ejA k f- g ki 

( 9ki — 9ki ) 



9k -i — 9ki — 9 hi 



9k -i - 9ki 



{Hk)i 



1 - {H k )i 



( 9ki - 9k% ) 



1 

1 — (H k )u 



( 9ki - fffct ) • 
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